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Tilting Saturn without tilting Jupiter: Constraints on giant planet migration 

R. Brasser^ and Man Hoi Lee^ 

ABSTRACT 

The migration and encounter histories of the giant planets in our Solar System can 
be constrained by the obliquities of Jupiter and Saturn. We have performed secular sim¬ 
ulations with imposed migration and A-body simulations with planetesimals to study 
the expected obliquity distribution of migrating planets with initial conditions resem¬ 
bling those of the smooth migration model, the resonant Nice model and two models 
with five giant planets initially in resonance (one compact and one loose configuration). 

For smooth migration, the secular spin-orbit resonance mechanism can tilt Saturn’s spin 
axis to the current obliquity if the product of the migration time scale and the orbital 
inclinations is sufficiently large (exceeding 30 Myr deg). For the resonant Nice model 
with imposed migration, it is difficult to reproduce today’s obliquity values, because the 
compactness of the initial system raises the frequency that tilts Saturn above the spin 
precession frequency of Jupiter, causing a Jupiter spin-orbit resonance crossing. Migra¬ 
tion time scales sufficiently long to tilt Saturn generally suffice to tilt Jupiter more than 
is observed. The full A-body simulations tell a somewhat different story, with Jupiter 
generally being tilted as often as Saturn, but on average having a higher obliquity. The 
main obstacle is the final orbital spacing of the giant planets, coupled with the tail of 
Neptune’s migration. The resonant Nice case is barely able to simultaneously reproduce 
the orbital and spin properties of the giant planets, with a probability ~ 0.15%. The 
loose five planet model is unable to match all our constraints (probability < 0.08%). 

The compact five planet model has the highest chance of matching the orbital and 
obliquity constraints simultaneously (probability ~ 0.3%). 

Subject headings: celestial mechanics - planets and satellites: dynamical evolution and 
stability - planets and satellites: formation 


1. Introduction 

In the last two decades, considerable progress has been made in regards to understanding the 
processes that shaped the dynamical evolution of the Solar System, and it seems the migration of 
the giant planets is a critical chapter in this story. Angular momentum exchange with a remnant 
planetesimal disc induces planetary migration, moving Jupiter slightly inward and Saturn, Uranus, 
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and Neptune outwards (|Fernandez &: id Il984l : iHahn &: Malhotral 119991 ) . The existence of many 
Kuiper belt objects in orbital mean-motion resonances (MMRs) with Neptune can be explained by 
the outward migration of Neptune. If the migration of the planets was smooth, the orbital eccen¬ 
tricities of the resonant Kuiper belt objects and the asymmetry in the resonant angle distribution 
of objects in the 2:1 resonance can be used to infer that Neptune has migrated by ~ 7AU on a 
timescale > 1 Myr and that the ini t ial mass of the planetesimal disc is ^ 50Mm, where M m is the 
mass of the Earth ( Malhotra 1995 : Hahn fc Malhotra 1999 : Murrav-Clav fc Chiang 2005 ). How¬ 
ever, this smooth migration model has difficulty reproducing several features of the Solar System, 
and alternative scenarios have been invoked. 


Thommes et al.l ( 19991 ) proposed a scenario in which Uranus and Neptune grew as cores between 


Jupiter and Saturn, were scattered close to their current locations due to planet-planet gravitational 
interactions, and had their random velocities damped by interactions with a remnant planetesimal 
disc. In this way the challenge of forming the ice giants at their current locations, with long orbital 
timescales and low (estimated) solid surface density at the time of formation, could be neatly 
side-stepped. Their work provided the first strong numerical evidence that the outer Solar System 
could have undergone considerable rearrangement from the time that the giant planets achieved 
their current masses. 

The successor of this model for the conf iguration of th e giant planets after the protoplane¬ 
tary gas disc has dissipated was presented by iTsieanis et al.l (j2005l j and is known as the Nice (or 
‘classic Nice’) model. It proposes that the giant planets were once in a considerably more com¬ 
pact configuration, in which Saturn started just inside the 2:1 MMR with Jupiter, and the ice 
giants were located at 11-13 and 13.5-17 AU. When the 2:1 resonance was crossed by the diver¬ 
gent migration of Jupiter and Saturn caused by the gravitational interactions with planetesimals, 
the system became dynamically unstable, exciting eccentricities and inclinations, and scattering 
Uranus and Neptune to close to their present locations (in some cases, after switching their order). 
This s cenario has many nice properties, including possibly explaining the Late H eavy Bombard¬ 
ment (jGomes et al.ll2nn5l ) and Jupiter’s Trojan population (jMorbidelli et al.l 120051 ). while typically 
preserving the regular satellite populations during the planetary encounters. 

However, even granting the classic Nice model its successes, there were obvious questions about 
the history: how did the system arrive at such an initial configuration? The timescales for migration 
induced by interactions with the protoplanetary gas disc to drive the giant planets i nto the inner 
system are shorter than the formation times for the object (e.g., iTanaka et al.ll2002l ). The recent 
revisions to the type I mig ration rate from a more sophisticated treatment of thermodynamics (e.g.. 


Paardekooper et al.ll2r)irilj have complicate d the situation, and it remains a challenge. Building upon 


work done by iMasset fc Sne Igrovd (1200111 and iMorbidelli fc Cridal (120071 ) , a variant of the classic 
Nice model was proposed by Morbidelli et al. ( 20071 ). This resonant Nice model attempts to bridge 
the gap between the gas-rich phase (typically studied with one or two planets embedded in a gas 
disc modelled using a hydrodynamics code) and the gas-free phase (studied with many planets and 
planetesimals using an Ai-body code). The authors use the observation that convergent migration in 
the gas-rich phase can result in capture into MMR and that a 3:2 MMR between Jupiter and Saturn 
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can prevent Jupiter from migrating inwards. They construct a number of stable configurations in 
which each giant planet is in resonance with its neighbours (with Saturn and the inner ice giant in 
the 3:2 MMR, and the ice giants in the 4:3 or 5:4 MMR). Subsequent migration due to interactions 
with the planetesimals breaks the resonances, and an instability resembling that of classic Nice 
occurs. 

Recent studies on giant planet migration have tried to dete rmine the actual evolu tion based on 
constraints from the secular arc hitecture of the gia nt planets (|Morbidelli et al.l 1200911. the curren t 
orbits of the terrestrial planets (jBrasser et al.l 120091 1 and the asteroid belt (jMorbidelli et al.ll2010l l. 
In essence, it is required that Jupiter scattered one of the ice giants outwards and thus migrated 
on a very short time scale to avoid both making the orbits of the terrestrial planets too eccentric 
and causing secular resonances to sweep through the astero i d bel t. Un fortunately th is scen ario 
had a very low probability of success, which led iNesvornvl ()201lh and iBatvgin et al.l (120121 1 to 
propose a scenario in which the Solar System began with five giant planets locked in MMRs. 


Nesvornv fc Morbidellil (j2012l l have studied a large number of hve planet systems to try to reproduce 


the current Solar System and match all the constraints, which they did with a probability of about 
5% in the best cases. 

One thing that has not yet been considered are the effects of migration on planetary obliqui¬ 
ties (the angles between the spin axes of the planets and their respective orbit normals) as the tilt 
angles of the giant planets preserve information about migration and encounter history. The case 
of Saturn is particularly interesting because of a coincidence between the spin precession rate of 
Saturn and the vertical secular eigenfrequency associated with Neptune, as well as the near match 
of Sat urn’s 27° obliquity t o tha t of Cassini state 2 (see Sectio n [2]) of the secular spin-orbit reso¬ 


nance. 


Ward &: HamiltonI (j2004l l and iHamilton &: WardI (j2004l l have proposed an eleg ant scenario. 


in w hich the secular spin-orbit resonance is resp onsible for tilting Saturii (but see iHelled et al 


20091 for a potential problem with this scenario). [Hamilton fc WardI (j2004l ) also note that, while 


they concentrate on changes in secular eigenfrequencies due to the changing gravitational potential 
of a gradually decaying Kuiper belt, all that matters is that the frequencies change and not the 
underlying mechanism, and so a change in frequency due to migration should also work. At the 
same time, even though Jupiter’s spin precessi on rate is close to the vertical secular eigenfrequency 
associated with Uranus (jWard &: Canupl l2006lb its 3° obliquity means that its obliquity evolution 
was very different from Saturn’s. 

Whether this spin-orbit resonance mechanism can work in the various migration scenarios is 
therefore a natural question. Here we perform secular and particle A^-body integrations and we 
initially consider smooth migration. Nice-like, and five planet models for the giant planets, and 
study the resulting obliquities. The obliquities of Jupiter and Saturn, in combination with the 
survival and orbital properties of four giant planets, provide strong constraints on the migration 
and encounter histories of the giant planets. In Section [2] we consider the conditions needed to tilt 
Saturn by spin-orbit resonance. In Section [3] we describe our numerical methods, and in Section 
0] our initial configurations. We present our results in Section [5l with the secular and iV-body 
simulations in Section 0 and Section [Ql The results are discussed in Section [6] and summarised 
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in Section [71 


2. Secular spin-orbit resonance and migration timescale 

We construct a fixed reference frame with unit vectors {i,j, k) where k is parallel to the total 
orbital angular momentum of the system (and hence perpendicular to the invariable plane) and i 
and j are in the invariable plane. Consider a planet whose unit vector of the spin direction s is 
tilted at an angle e from the unit orbit normal n, i.e., cose = s ■ n. The orbit has an inclination 
cosi = n ■ k. If there are no perturbations on the orbit and thus n is fixed in space, the torque 
from the Sun on the oblate figure of the planet (and on any satellites whose orbits are locked to the 
planet’s equator plane) causes s to process uniformly about n at the frequency q ;(1 — cose. 

Here e is the orbital eccentricity and a is the precession constant, given by 

^ 3GMq ( J 2 + q' \ 

2i'a^ \ C + i / ’ 


where G is the gravitational constant, Mq is the solar mass, a is the orbital semi-major axis of the 
planet, i' is its rotation frequency, J 2 is the second zonal harmonic coefficient of its gravity field, C 
is its moment of inertia normalised to mB?, m and R are the mass and radius of the planet, £ is 
the orbital angular momentum (normalised to mR^v) of the satellites whose orbits are locked to 
the planet’s equat or, and q'IJ 2 is the ratio of the torque on the satellites to that directly exerted 
on the planet (see I Ward fc Hamiltonll2nn4l for explicit definitions of £ and q'). 

Secular perturbations from other planets force n to process about k. In the simple case where 
n processes at a constant rate s, there are either two or four locations where s remains coplanar 
with n and k, and s and n process at the s ame rate s about k. These secular spin-orbit resonance 


conf igurations are calle d the C assin i states ([Colombo 


1966; 


Pealelll969l . [l97J). 


Ward fc HamiltonI (j200J) and [Hamilton fc WardI (|2004l i suggest that Saturn’s obliquity of 


nearly 27° may be the result of Saturn being caught in Cassini state 2 with the vertical secu¬ 
lar eigenfrequency sg associated with Neptune0 They show that Saturn’s obliquity eg can be 
increased from an initially small value to 27° by capture into Cassini state 2 if [sg/agl > 1 initially 
and decreases to < 1. They also show that the more complicated nodal regression of Saturn caused 
by perturbations from Jupiter and Uranus do not affect the capture into spin-orbit resonance and 
the associated growth in Saturn’s obliquity. 

One way to decrease [sgl (and hence jsg/asj) is through the outward migration of Neptune. 
The time scale r* = ag/sg for the decrease in |sg| must be adequately long for Saturn’s spin to 
stay in resonance. The critical time scale Ts^min can be estimated from requiring the time that 
it takes for the obliquity of Cassini state 2 to move by an angle comparable to the width of the 
separatrix to be longer than the libration period about Cassini state 2. With the obliquity of 


^We follow Bretagnon 1 1974 ') and Laskar |[l990l) an d den ote t he vertical secular eigenfr equencies associated with 
the Rian t planets by se, S 7 , and ss- IWard fc HamiltonI ([2004 )_aad Hamilton_&Wajd (l2004l ~) use the notation gie, gn, 


and gis’, Applegate et al. ( lOSfill 


use ge, gr 


and gs’, while Murray fc Dermottl 1 199 £ ) use fe, R, and /g. 
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Fig. 1.— (a) Evolution of Saturn’s obliquity e in a set of simulations in which Saturn’s orbit normal 
n is inclined by i = 0.1° from the normal to the invariable plane k and forced to process about 
k at rate s. The orbital eccentricity e = 10“^, the spin precession constant as = 0.79"yr“^, and 
s changes linearly from —5as to —0.88as according to s{t) = (—5 + t/Tg)as- Saturn remains in 
Cassini state 2 to the end of the simulation if Ts > ~ 700 Myr or escapes from the resonance 

before the end of the simulation if Ts < Ts^min- (b) Critical time scale rs,min as a function of orbital 
inclination i for both Saturn and Jupiter (with aj = 2.67" yr“^). The dashed lines show Ts^min 
according to Equation ([2]) with the final obliquity e ~ 29°. 


Cassini state 2 given by coses ~ — (s8/«s)cosi (when es is large), the half width of the separa- 
trix Aes = 2(tani/tanes)^/^, and the libr ation frequency cjiih = (—asss sines sinfor small- 
amplitude libration abo ut Cassini state 2 (jHenrard &: Murigandd 119871 : IWard k. HamiltonI l2004l : 
Hamilton k, Ward 2004 1. Then 2Aes/es ~ 27r/a;iib gives 


vr 


2as sin i sin es 


( 2 ) 


The condition ss ~ fomid bv [Hamilton fc WardI (120041 ) gives a similar Ts^min (but without the 
coefficient vr/2). Boue et al. ( 2009l l have also shown using a different analytic argument that the 
critical time scale Ts^min oc l/(asi). 

We have verified Equation m with numerical simulations using t he symplectic i ntegrator 
SyMBA, which was modified to include the evolution of the spi n axis (ILee et al.l 120071 : see also 
Section 3). We integrated Saturn alone with as = 0.79" yr“^ ( Tremain^ 1991 1. e = 10 ^ and 
imposed orbital nodal precession rate s that changes linearly from —5as to —0.88as according to 
s{t) = (—5 + t/Ts)as- The spin axis is initially parallel to k. Fig. [TKa) shows the evolution of 
Saturn’s obliquity in a set of simulations with i = 0.1° and different r^. Saturn is captured into 
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Cassini state 2 and remains there until the end of the simulation (when es is slightly larger than 
27°) if Ts > Ts^min ~ 700 Myr. For Tg < Ts^min) Saturn is initially captured into Cassini state 2, but it 
escapes from the resonance when the obliquity becomes too large; the final obliquity increases with 
increasing Tg and can be more than 10° for Tg > From similar simulations with i = 0.01°- 

10°, we have determined Ts^min as a function of i, and the results are shown in Fig. HKb) for Saturn, 
as well as Jupiter with aj = 2.67" yr“^. We confirm Equation ([2]) over a wide range in i and find 
that for Saturn ~ 20 Myr for z of a few degrees. We note that 20 Myr is roughly the upper 

l imit of the typical m igration time scale of the giant planets caused by planet esimal scatterin g 
(|Morbidelli et al.ll2010l b Also, in our earlier simulations of the classic Nice model (|Lee et al.ll2007l b 
we found that the orbital inclination of Saturn can reach ~ 5° (and those of the ice giants 10-15°) 
during the encounter phase in some simulations (see, e.g., Fig. 2 of iLee et al.ll2007l ). 

Summarising, it appears possible to tilt Saturn to its current obliquity if Iss/asI changes 
slowly enough. In the next sections we investigate whether or not this can be realised in the 
various migration scenarios. 


3. Numerical methods 


To study the obliquity evolution in the various migration scenarios, we performed several types 
of numerical simulations of the Sun and the giant planets. One set of simulations consisted of inte¬ 
grating the secular evolution of the migrating giant planets without encounters and mean-motion 
resonances. The second set of simulations integrated the full Wbody proble m with planetesimals . 
These integrations used the Kepler-adapt ed symplectic Wbody cod e Sv MBA (IDuncan et al.lll998l ') , 
a descenda nt of the or i ginal techniques of lWisdom fc HolmanI (|l99ll ) and iKinoshita et al.l (|l99ll ). as 
modified bv iLee et al.l (120071 ) to incorporate spin evolution. There are some additional modifications 
for both the secular and Wbody simulations that we now describe. 


3.1. Secular Simulations 


Since we expect the obliquity evolution to be driven by resonances between the nodal re¬ 
gression and spin precession frequencies, both rather small compared to the orbital frequencies, a 
useful approximation to the full Wbody behaviour for exploratory simulations can be obtained by 
modelling the p urely secular evolution of the system. We have adapted the spin-SyMBA code of 
Lee et al.l (120071 ) to remove mean-motion effects and close encounters when desired. The method 
is the following. 

We consider a system made up of the Sun and some number of planets. The planets are 
assumed to undergo principle axis rotation, and the effects of the oblateness of the planets on the 
orbital evolution are neglected. The integration step has the form 

^T/2^r/2^T/2^r^r/2^r/2^r/2 (3^ 


where r is the time step and each term is an operator. The operator D advances the planets along 
their osculating Kepler orbits, S advances the spin axes according to the mutual torques between 
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planets and the torques due to the Sun, R handles the secular interactions between the planets, 
and M generates radial m igration (whose r ecipe is discussed later in this subsection). Unlike with 
the full code described bv iLee et al.l (|2007l b the integration of the secular system is performed in 
heliocentric coordinates rather than canonical helioc entric, and so there is no linear drift term to 
correct for the barycentric momentum of the Sun ([Duncan et ahl Il998l i . Since we do not treat 
close encounters in the secular approximation, we do not need the recursive subdivision of the time 
step of SyMBA. The Kepler drift operator is straightforward, and the spin torque operator follows 
Lee et ahl (|2007l L but the secular operator replacing the interaction term is new. 

The secular interaction between the planets can be treated in several ways. Here we forgo 
the eigenvalue approach and instead directly integrate Lagrang e’s planetary equations with the 


first-order secular potential (see, e.g., iMurrav &: DermottI Il999li . Let a, e, i be the semi-major 


axis, eccentricity, and inclination of a planet, zu the longitude of perihelion, H the longitude of the 
ascending node, and n the mean motion. To apply the R operator, for each pair of planets (with 
indices j, k; j 7 ^ k) we write the disturbing function on planet j due to k as 


Rj}^ Ojj 


1 


9 COs(tx7j -(- ^jj U U COs(Hj ^k) 


(4) 


Here the matrix elements A and B are functions of the m asses and semi-major axes and they 
contain information about the coupling between the planets ([Brouwer &: van Woerkomlll95[ll ). We 
can then compute the time derivatives 


1 dR 


^jk — 
VJjk 


jk 


1 OR 


rij ttj Cj 


^ “ i dwi 


dR 


■jk 


njajej dej 


I'jk — 


il,jk 


■jk 


Tij aj ij dflj 


1 dR 


'jk 


rij a - ij 


dii 


(5) 

( 6 ) 


The total rate of change in eccentricity of planet j is given by 


e-j — ^ Cjk 


(7) 


and similarly for the changes in the other elements ij , ujj , and . These equations are integrated 
using a fourth-order Runge-Kutta scheme. 

It may seem odd to evolve what is fundamentally a secular ring problem using point particles, 
as this requires repeatedly solving Kepler’s equation during the Kepler drifts, which should be 
unn ecessary in a sec ular model. However, in this way the well-tested spin vector infrastructure of 
the iLee et ahl (|2007l j modifications to the SyMBA code could be reused. We have confirmed by 
comparing with a different secular code (one which solves the corresponding eigenproblem instead) 
and with the full A^-body integration that this approach recovers the expected evolution. 

Note that the secular step itself is rather time-consuming compared to a simple force evaluation 
due to the need to compute multiple trigonometric functions and the Laplace coefficients. The 
separation between the orbital time scale and the secular, spin, and migration time scales provides 


























several opportunities for optimisation. Advancing the secular term R every (orbital time scale) 
step is unnecessary, and applying it every 10 steps serves to make the runtime comparable to the 
A^-body run, without changing the resulting evolution. Another improvement can be obtained by 
dividing each of the secular, spin, and migration time scales by some constant factor (making sure 
to keep them several orders of magnitude longer than the orbital time scale) during the simulation 
and then scaling the resulting output. This is possible because the behaviour of the spin-orbit 
resonance depends not on the absolute values of these time scales but on their ratios. In our 
regime, using a factor of 10-20 produces negligible errors relative to the original simulation and 
provides almost all of the possible speed benefits. We have verified that our results are insensitive 
to these approximations, which are used only in the secular integrations. 

In the standard picture, the final configuration of the giant planets should evolve mostly due 
to mutual interactions and interactions with an external remnant planetesimal disc (such as the 
Kuiper belt). Accordingly, with the possible exception of objects which get scattered out into the 
Kuiper belt and Scattered Disk (jPuncan &: LevisonI 119971 ). imposing a prescription for migration 
and random velocity damping during the secular simulations should be a reasonable approximation 
to the effect of the planetesimal disc. Let t be th e time since the start of the simulation, and r the 
migration time scale. We follow iLee et al.l (120071 ). and take as migration expressions 


/ j. / \ 

Ofc =-exp(-t/r). 


r 


and 


a-k = 


Attk {t 


exp[-tV(2T^)], 


( 8 ) 


(9) 


where Aa^ is a parameter measuring the amount of distance the planet should travel due to the 
migration to end up on its current orbit. In the secular problem, Aa^ is simply the difference 
between the current value and the starting value. Equation ([8|) results in faster migration at earlier 
times than Equation Q, but the two resulting trajectories meet again at t = 2r when the distance 
travelled is ~0.86Aafc, after which Equation ([9]) has faster migration. Equation ([8]) has a continually 
decreasing d, whereas Equation ([9j) has d increasing until a maximum at r (corresponding to a 
maximum speed of ~ O.GAofc/r) beyond which it decays to zero. The M operator of Equation 
([3]) simply applies the changes in a to the planets. As in iLee et al.l (j2007li . we evolve all planets 
under Equation ([8|) in our integrations of the smooth migration model, and we use Equation ([HI) for 
Jupiter and Saturn and Equation ([9|) for Uranus and Neptune in our integrations of the Nice model. 
Note that we do not apply any eccentricity or inclination damping to the secular simulations: as 
there is no stirring mechanism, e and i would damp too quickly to serve as a useful model for the 
spin-orbit resonance crossing. We evolve every simulation to lOr. 

This concludes our description of the secular simulations. We now move on to discuss the full 
A-body simulations. 
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3.2. AT-body simulations 

In the full model the planets migrate because of their interaction with the planetesimal disc 
and by mutual encounters. We integrated the giant planets and a planetesimal di sc consisting of 
several thousand planetesimals using the spin-SyMBA code from iLee et al.l (|2007l b We made an 
additional modification to the code so that it only outputs the orbital elements of the planets and 
not those of the planetesimals in the disc, and stops the integration once the number of planets is 
fewer than four. The i nitial conditions for the s e sim ulations are described in the next section and 
are similar to those of iNesvornv &: Morbidellil (|2012l i. The outer Solar System and the obliquities 
of the giant planets were integrated with a time step of 0.35 yr. Planets and planetesimals were 
removed once they were further than 1000 AU from the Sun (whether bound or unbound) or when 
they collided with a planet or ventured closer than 0.3 AU to the Sun. We evolve every simulation 
to 500 Myr or until fewer than four planets remain, whichever occurs first. The simulations were 
run on the TIARA computer cluster at the Institute for Astronomy and Astrophysics, Academia 
Sinica, in Taiwan. 

For these A-bocW simulations, a minor adjustment to the SyMBA algorithm as implemented 
in the Swift packagqj is needed. SyMBA decomposes the Hamiltonian in canonical heliocentric 
coordinates i.e. heliocentric positions and barycentric velocities, and employs a multiple time step 
technique to handle close encounters. When we first integrated our initial conditions for the resonant 
Nice configurations in the absence of planetesimals, we found unexpectedly that no more than half 
of the initial conditions were stable for 1 Gyr when we used SyMBA, whereas they were all stable 
under both the original Wisdom-Holman method and a pure canonical heliocentric integrator i.e. 
without any treatment of close encounters. This instability is caused by Jupiter and Saturn being 
in the 3:2 resonance, where they are close enough at conjunctions that they are considered to be 
undergoing encounters by the SyMBA algorithm, and their time steps are subdivided at every 
conjunction. The time step subdivision is achieved in the SyMBA algorithm by using a partition 
function to decompose the gravitational force between two planets into forces that are non-zero 
only between two cut-off radii, Rk+i < r < R^. The simplest partition function is the {2^ -\- l)th 
order polynomi a.! in x that has //(O) = 1, /^(l) = 0, and all derivatives up to the Hh derivatives zero 
at X = 0 and 1. iDuncan et al.l (| 19981 ) found that the third-order polynomial fi{x) = 2x^ — 3x^ -|- 1 
worked well for many situations, which did not include repeated encounters on orbital time scales 
over 100 Myr-1 Gyr. For this more challenging situation, a smoother partition function is needed, 
and the next appropriate polynomial is / 3 (x) = 20x^ — 70x® -|- 84x^ — 35x^ -|- 1. We found that 
the use of fz{x) sufficed to keep each of the resonant Nice initial conditions stable under SyMBA 
for 1 Gyr. Of course, in our actual simulations with planetesimals, the planets will have migrated 
on time scales much shorter than the numerical instability time scales, but it is useful for testing 
purposes that the integrator does not fail on the time scales of interest. 


^Available at http://www.boulder.swri.edu/'hal/swift.html 
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4. Models and initial conditions 


We restrict o urselves to four classes of initial conditions, namely a version of the smooth mi- 
gration scenario oflHahn &: Malhotral (|1999l l with four planets, a variant of t he resonant Nice model 
( Morbidelli et al. 2007l i wi th four planets, and two f ive planet models from Nesvornv &: Morbidelli 
(|2012l b We integrate the iHahn &: Malhotral (|l999l l initial conditions with the secular equations 
only, the resonant Nice with the secular equations and full A^-body model, and the five planet cases 
with full Wbody only. 


4.1. Smooth Migration 

We set the initial semi-major axes of Jupiter and Saturn to be aj = 5.40 AU and as = 8.73 AU, 
respectively. Uranus and Neptune are given random semi-major axes uniformly drawn from the 
ranges 14.19-18.19 AU and 21.07-25.07AU. This configuration means that Saturn begins just 
beyond the 2:1 resonance with Jupiter at 8.6 AU. Eccentricities are set to 0.001. 


4.2. Resonant Nice Confignration 

In this case, the four giant planets are all initially in resonance. Saturn and Jupiter are in the 
3:2, as are Saturn and the inner ice giant. Th ere are two distinct stable configurations for Uranus 
and Neptune found bv iMorbidelli et al.l (j2007l i. namely the 4:3 and the 5:4; we study only the 4:3, 
which has the out er ice giant at 11 . 6 AU . To generate the initial conditions, we used an approach 
similar to that of IMorbidelli et al.l (|2007l h but instead of a hydrodynamic code, we employed the 
A^-body code SyMBA, with forced migration and eccentricity damping to model the effects of the 
proto-planetary gas disc. We began with Jupiter and Saturn and placed them just outside the 
3:2 mean-motion resonance. Jupiter was forced to migrate outwards and Saturn inwards, and the 
migration rates were chosen so that they would be nea rly stationary a fter capture into the 3:2 
resonance. The eccentricity damping rate is specified by ()Lee et al.l 120071 ) 




( 10 ) 


We adopted Kp = 100 for Jup i ter an d adjusted Kp for Saturn to match the equilibrium eccentricities 


mi 

reported bv IMorbidelli et al.l (j2007l i. After the Jupiter-Saturn pair had reached equilibrium in the 
3:2 resonance, we added an ice giant just outside the 3:2 resonance with Saturn and forced it to 
migrate inwards until Saturn and the ice giant reached equilibrium in the 3:2 resonance. This 
process was repeated for the 4:3 resonance between the ice giants. For the ice giants, which should 
be undergoing type I migration, we assumed that the values o f Kp a re identical and adjusted Kp 
until we matched the eccentricities reported bv iMorbidelli et al.l (120071 ) at each stage of this process. 
To model the dispersal of the gas disc, we continued the integration with an exponential decay in 
the migration and eccentricity damping rates. 
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4.3. Five-Planet Configurations 

For the five planet cases we have chose n two configurations which br acket the most compact 
and least compact configurations studied bv iNesvornv &: Morbidellil (j2012l i. The configurations we 
studied were 3:2, 3:2, 4:3, 4:3 and 3:2, 3:2, 2:1, 3:2. In the first case the outermost ice giant was at 
14.2 AU while in the latter case it was at 22.2 AU. The initial conditions for the first (compact) 
configuration were generated in a similar manner to the resonant Nice configuration, but with the 
addition of a third ice giant in 4:3 resonance with the second one. The initial conditions for the 
second (loose) configur ation were provided by D. Nesvorny. Th ese two configurations br acket the 
resonant Nice model of iMorbidelli et al.l (|2007l i and the model of lHahn &: Malhotral (|l999l i in terms 
of spacing between the planets, but with an additional ice giant. 

For the initial conditions of the A^-body simulations with planetesimals, we selected randomly 
from the last five frames of output from the gas dispersal integration for the resonant Nice and 
compact five planet cases, and the system is rescaled so that initially Jupiter’s semi-major axis 
aj = 5.45 AU. We confirmed that the last five frames are stable for 1 Gyr (in the absence of forced 
migration and damping), after we made the reported minor adjustment to the SyMBA algorithm 
(see Section [3.21) . For the initial conditions provided by D. Nesvorny we generated five frames by 
integrating the planetary system without the planetesimal disc for 500 yr and output the orbital 
elements every 100 yr. 

The planetesimal disc consisted of 2000 planetesimals for the resonant Nice and compact five 
planet cases and 3000 for the loose five planet case. The surface density of the planetesimal disc 
scales with heliocentric distance as S oc r~^. The outer edge of the disc was always set at 30 AU. 
The inner edge was A = 0.5 AU from the outermost ice giant in the resonant Nice and compact 
five planet cases and A = 1 AU for the loose five planet case. The total disc mass was 5OM0 for 
the r esonant Nice model, 35114^ for the compact five planet case and 20M^ for the loose five planet 


case. iNesvornv &: Morbidellil (|2012li have found that a disc of ~ 5OM0 gives the best results for the 
specific resonant Nice configuration that we adopted. We decreased the disc mass for the compact 
five planet case to account for the mass of the extra ice giant. All configurations have roughly the 
same surface density in planetesimals. We generated 256 different realisations of each frame by 
giving a random deviation of 10“® AU to the semi-major axis of each planetesimal, for a total of 
1 280 simulations per system. 

In all cases, we u sed the cu r rent m asses and spin precession constants of the giant planets, with 
the latter taken from iTremaine ( 1991 ). For the current orbital semi-major axes, the spin precession 
constants are aj = 2.67"yr“^, as = 0.793"yr“^, au = 0.0448"yr“^, and on = 0.0117"yr“^ for 
Jupiter, Saturn, Uranus, and Neptune, respectively. There is an ambiguity about interpreting the 
ice giants’ values in the resonant Nice and five planet cases, because the ice giants can switch 
places and, in the five planet cases, any one of the ice giants can be ejected. However, the masses 
of Uranus and Neptune are very similar (14.5M0 and I 7 .IM 0 ), and in an odd coincidence, their 
spin precession constants are also very similar when scaled to a common semi-major axis, so the 
ambiguity presents much less difficulty in practice than might have been expected. Thus we adopted 
the parameters of Uranus and Neptune for the initially inner and outer ice giants in the resonant 
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Nice model, and added a third ice giant with a mass of 15.(the average mass of Uranus and 
Neptune) and the precession constant of Neptune for the compact five planet configuration. For the 
loose five planet configuration provided by D. Nesvorny, the masses of the initially inner, middle, 
and outer ice giants were equal to ISM^, the mass of Uranus, and the mass of Neptune, respectively, 
while the precession constants were the same as in the compact case. We set the initial spin vectors 
of the planets to be only slightly perturbed (0.001 radians) from the inertial z-axis. The initial 
orbital inclinations are specified below for the different types of simulations. All other angles are 
set at random values. 


5. Results 

In this section we present the results of our secular and A^-body simulations. We discuss the 
smooth migration case first to explain the underlying mechanism of how to tilt the gas giants, 
followed by the resonant Nice and five planet simulations. 


5.1. Secular Simulations 

To understand how the secular spin-orbit resonance scenario works in the full planetary mi¬ 
gration context, but without the complicating factors present in noisy A^-body simulations, we 
apply the secular spin approach of Section IS. II as an intermediate scheme between the one-planet 
simulations in Section [2] and the full A^-body runs to two of the histories (smooth migration and 
resonant Nice) described in Section [H 


5.1.1. Smooth migration 

The simplest model to consider is a purely secular model of the smooth migration scenario, as 
it behaves as the analytic theory would predict. Figure [2] presents an overview of a successful run 
that reproduces the current obliquities of the gas giants. In this example r = 32 Myr and initial 
i = 2° for all planets. The upper panels describe the evolution of the orbits, with the migration 
in semi-major axes imposed by Eq. ([8]). The inclinations oscillate around the forced values due to 
secular interactions. The lower panels describe the evolution of the system spins. The obliquities e 
of all four planets remain low until approximately 20 Myr, during which time the evolution of e is 
driven entirely by changes in the inclination, as can be seen by comparing the bottom left and top 
right panels. However, after about 20 Myr the obliquity of Saturn rises to a mean of ~ 33° with a 
superimposed oscillation of about 4°. That this obliquity rise is due to capture into a secular spin- 
orbit resonance is clear from the bottom right panel, which compares the spin precession frequencies 
a of the planets with the three non-degenerate vertical secular eigenfrequencies |s6|, [syl, and |s8| 
of the system. At the beginning of the simulation, no eigenfrequency is particularly close to any 
spin precession frequency of a planet, with the slowest eigenfrequency [sgl (which is that dominated 
by Neptune’s contribution) situated between the spin precession frequencies of Jupiter and Saturn. 
However, as the semi-major axes change due to migration, so do both the eigenfrequencies and the 
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Fig. 2.— Evolution of a successful purely secular smooth migration model with migration timescale 
T = 32Myr and initial inclinations i = 2°. The upper left panel shows the orbital semimajor axis 
a, perihelion distance q, and aphelion distance Q for Jupiter (red), Saturn (green), Uranus (blue), 
and Neptune (magenta). (Note that here the eccentricities are low enough that the lines cannot 
be distinguished.) The upper right and lower left panels show the orbital inclinations i and the 
obliquities e, respectively. In the lower right panel, the coloured lines correspond to the spin 
precession frequencies a of the planets, and the three black lines are the nondegenerate vertical 
secular eigenfrequencies [sel (highest), |s 7 | (middle), and jssl (lowest). Note that |s8| is initially 
between the spin precession frequencies of Jupiter and Saturn. 
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Fig. 3.— Resulting obliquities for Jupiter (ej) and Saturn (es) in secular simulations of the smooth 
migration model as a function of r. Each colour corresponds to a different initial inclination, from 
0.01° on the left to 2° on the right. The diagonal cross marks the observed configuration of Jupiter 
and Saturn. 
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spin precession frequencies. After 20Myr, when Saturn is observed to begin tilting over, we see 
that Isgl matches the spin precession frequency of Saturn, and it is this secular spin-orbit resonance 
which could explain the current tilt. 

Not every set of parameters results in such a Saturnian obliquity: as discussed in Section [2l 
both r and i play a role such that r x i must exceed a certain value. Fig. [3] shows the resulting 
Jovian and Saturnian obliquities from simulations with r = 1-100 Myr and initial i = 0.01-2°. In 
general, as the time scale increases so does the resulting Saturnian obliquity, up to a maximum of 
~30°. This is very close to the observed 27°, thus providing a satisfying explanation for Saturn’s 
obliquity. Moreover, in this migration time scale regime, at least an inclination of 1-2° is required, 
which also gives a Jovian obliquity compatible with the current value. Thus smooth migration can 
work at large r x i (> 30 Myr deg): the minimum r x z producing a simulation with hnal Saturnian 
obliquity > 15° is 32 Myr deg, although some simulations with r x z above this limit have final 
Saturnian obliquity < 15°. This joint constraint will prove to be useful. 

Although the smooth migration model is able to reproduce Saturn’s obliquity if r x z is suffi¬ 
ciently large, we want to emphas ise that it has great d ifficulties in explaining the c urrent secular 


archi tecture of the giant planets (|Morbidelli et al.ll2009l ) and the terr estrial planets (IBrasser et al 


20091 ) ■ as well as the orbital properties of the main belt asteroids (jMorbidelli et al 


201^ and 


should not be considered as a proxy for the past evolution of the Solar System. 


5.1.2. Resonant Nice 

In principle the situation with the resonant Nice model should be very similar. The outward 
migration of Saturn and the ice giants will lower the spin precession frequency of Saturn slightly and 
the eigenfrequencies considerably, so the same resonance crossing should occur, and it does. The 
difficulty becomes apparent when we examine an equivalent to the smooth migration run shown in 
Fig. [21 

Fig. m shows the evolution of a resonant Nice model with r = 10 Myr and initial zn = 4° 
(the initial inclinations of the other planets are zero). This relatively high inclination of Neptune 
was chosen to mimic the typical inclination value during the encounter phase from the A^-body 
simulations. The Nice initial configuration is significantly more compact, which has the effect of 
increasing the eigenfrequencies relative to the spin precession frequencies. In particular, the same 
frequency that successfully tilted Saturn in the smooth migration scenario, [sgl, is now larger than 
Jupiter’s spin precession frequency and must cross it to reach its final value near Saturn’s. The 
situation for Jupiter is exactly comparable to the case of tilting Saturn, and the only question is 
whether the product of the migration time scale and the inclination is large enough at the time of 
the resonance crossing near 10 Myr to produce a significant Jovian obliquity. (Recall that Jupiter 
can be tilted by a faster decrease in [ssl than Saturn; see Fig. [TJ) Unfortunately, the answer is 
yes. In this run, at the end we have ej ~ 6°, which is much greater than the observed 3°, even 
though Jupiter does not stay in resonance until the end. Moreover, Saturn’s obliquity of ~ 8° is 
not substantially larger. To increase es we require an increase in either zn or r or both, but since 
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Fig. 4.— Same as Fig. [2l but for a secular resonant Nice model with r = lOMyr and = 4°. 
Note that jssl now crosses the spin precession frequency of Jupiter before it reaches that of Saturn, 
resulting in Jupiter tilting to ~ 6°. 
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Fig. 5.— Same as Figj3l but for secular simulations of the resonant Nice model with zn = 1-10°. 
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the same resonance is at work and the crossing times are only separated by ~ r, improving the 
Saturn value will make Jupiter’s obliquity unrealistically large. Fig. [5] summarises the outcomes 
from a series of resonant Nice runs with r varying from 1 to lOMyr and the initial zn from 1-10°. 
As feared, ej ~ es- the mean and standard deviation of ej/es = 0.67 ± 0.12, whereas the observed 
ratio is 0.11. The only runs which show eg > 15° have 10° < ej < 20°, and have long migration 
times (8 or lOMyr). 

However, the secular models leave out many relevant dynamical effects. For example, we apply 
semi-major axis migration via our prescriptions, and the real system will undergo a more sudden 
transition during resonance crossing or breaking and will be far more stochastic. Perhaps Jupiter’s 
spin-orbit resonance crossing happens much more quickly than Saturn’s, preventing too large a 
tilt. In the secular simulations we also let the inclinations undergo the usual evolution. Some 
experiments (not presented here) introducing inclination damping proved unfruitful as the Jupiter 
crossing occurs first, and lowering the inclination enough to avoid tilting Jupiter results in a far 
too small Saturnian obliquity. In the real system, one could imagine that the inclinations are low 
during the Jupiter spin-orbit resonance crossing and higher due to encounters when Saturn crosses, 
recovering the smooth migration model behaviour. To determine whether these possibilities are 
realistic scenarios we must move beyond the secular models to more realistic N-body models, and 
to those we now turn. 


5.2. Al-body Simulations 

We saw from the secular simulations that in the resonant Nice model the obliquities of Jupiter 
and Saturn are similar because of Isgl crossing the spin precession frequencies of both planets. 
This did not occur in the smooth migration setup because of the wider spacing of the planets|§ 
We now present the results of our full iV-body simulations. In what follows we need to quantify 
what we consider to be a successful simulation that adequately reproduces the secular properties 
and spacing of the outer Solar System, in addition to the obliquities of Jupiter and Saturn. To 
that end we id entify several conditions that the system must adhere to. These are similar to the 
requirements of iNesvornv &: Morbidellil (j2012l ') but with some modihcations. First, we need to have 
four giant planets (Jupiter, S aturn , and two ice giants) at the end of the simulation - criterion 
A of iNesvornv &: Morbidellil (j2012l '). Second, the planets need to be within 10% of their current 
orbits, i.e., 0.9 < afjac < 1.1 where af is the hnal semi-major axis of each planet at the end of the 
simula tion, and Uc is its current semi-major axis - similar to criterion B of iNesvornv fc Morbidelli 
(j2012l ) which imposes a margin of 20%. Third, we impose a restriction on the difference in the 
longitudes of perihelion, Arojg = zuj — rug. The current secular evolution of the eccentricities 
of Jupiter and Saturn is such that Arojg circulates throug h 360° and that t he ec centricity of 
Jupiter (Saturn) reaches a maximum when Arojg = 180° (0°). iMorbidelli et al.l (j2009l ) have shown 


®Although the classic Nice initial configurations ( Tsiganis et al. 2005ll are also less compact than the resonant 
ones, they have |s8| larger than Jupiter’s spin precession frequency and have the same problem in the secular models. 
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that, at the very least, encounters between an ice giant and Saturn are needed to produce this 
secular behaviour. Thus, we require that Avjjs circulates. We used Fourier analysis over the last 
5 Myr of the simulation to es tablish the behaviour of Arojg. This approach differs slightly from 


Nesvornv &: Morbidellil (j2012l l who keep track of the amplitude of the eccentricity eigenmode of 


Jupiter (their criterion C). For the final obliquities, we require that £j < 5° and eg > 15°. 

During the course of this investigation, we became aware of Jupiter often being too close to 
the Sun at the end of our Wbody simulations. For example, the mean and standard deviation of 
the final aj = 4.98 ± 0.21 AU for the simulations of the resonant Nice model with four planets at 
the end. This affects not only whether the planets are within 10% of their current orbits, but also 
potentially the obliquities of Jupiter and Saturn, because the secular spin-orbit resonance problem 

is not scale free. The precession constants scale with the semi-major axis of Jupiter as a oc 

—3/2 

while the secular eigenfrequencies scale as s oc Oj , if we fix the semi-major axis ratios. So the 
final \s/a\ would be smaller by ~ (3/2)Aoj/oj if the final aj is closer than the actual value 5.20 AU 
by Aaj. Since the initial aj is in some sense arbitrary, we rescale each simulation so that the 
final aj = 5.20AU. If the final Isg/agI < 1 (i.e., there is resonance crossing) before rescaling and 
Iss/asl > 1 (he., there is no resonance crossing) after rescaling, we use the obliquity of Saturn 
before the crossing in the original (unsealed) simulation as its final value. We apply the same 
procedure to [sy/ajl for Jupiter. 


5.2.1. Final Orbits and Obliquities 

First we analyse the outcome from the A-body simulations of the resonant Nice model, with 
the planets residing initially in the 3:2, 3:2, 4:3 resonant chain and a planetesimal disc of 50M^. 
In Fig. [6] we plot the final ej versus eg. The values of both ej and eg are averages taken over the 
last 50 Myr in order to wash out any large-amplitude, long-period oscillations, except, as noted 
earlier, when there is secular spin-orbit resonance crossing before rescaling but no crossing after 
rescaling (|s 7 /aj| for Jupiter and jsg/agl for Saturn), in which case the obliquity is taken before 
the crossing. The left panel is a log-log version of the right panel. The small black dots depict the 
final obliquities for simulations that have four planets (but the planets may or may not be near 
their correct positions). The larger blue dots refer to simulations where the planets are also near 
the correct positions (but Arujg may librate or circulate). Finally, the red dots denote simulations 
that match all of our criteria in orbital properties for success. The dotted line denotes ej = eg. 

There are several visible trends in Figure [H First, the obliquities of both Saturn and Jupiter 
take on values ranging from under 1° to over 30°. They generally do not follow the line ej = eg. 
Second, there are roughly equal numbers of cases with ej > eg and ej < eg. We find that ej > eg 
50% of the time for all simulations that produce four planets (irrespective of their final orbits), and 
ej > eg 51% of the time for the cases when the planets are close to their current orbits. However, 
the mean and standard deviation of ej/eg = 1.86 ±2.75 and 2.24±3.40 for each sample, indicating 
that the distribution of ej/eg has a longer tail at large values. Compared to the secular simulations 
of the same model in Section [5. 1.21 the mean ej/eg is even further from the observed value of 0.11, 
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Fig. 6.— Final obliquities of Jupiter and Saturn in the A^-body simulations of the resonant Nice 
model. Each panel depicts the same outcome but with a different scale. Black dots denote simula¬ 
tions with four giant planets. Blue dots denote outcomes with four planets in the right positions. 
Finally, red dots denote outcomes that match all our orbital criteria. The asterisk marks the 
observed configuration of Jupiter and Saturn. 


but there is a much larger spread in ej /es • 

Unfortunately only very few simulations matched our criteria for being considered successful. 
A total of 84 runs out of 1280 (6.6%) had the planets in the right positions. Of these only 13 
have Arojs circulate. Of these 13 only one has es > 15° and ej < 5° (see also Tabled!). However, 
i t is n ot clear how essential of a constraint the circulation of Atujs is. iNesvornv fc Morbidelli 
(j2ni2l ) state that requiring the circulation of Atujs is the real bottleneck in matching the current 
orbital architecture. Removing this restriction yields 2 additional good cases matching the planetary 
architecture and the obliquities. In summary, it appears that the 4-planet model has some difficulty 
reprod ucing the current Solar System. This fact was already established bv INesvornv fc Morbidelli 
(j2012l l. who reported similar statistics for the orbital constraints, but the probability is even lower 
(~ 0.08-0.23%) when we consider the obliquities. 

Next we report on the outcome of the A-body simulations with five initial planets. Figure d 
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Fig. 7.— Same as Fig. [H but for the A^-body simulations of the compact (3:2 3:2 4:3 4:3) five 
planet model. 

is similar to Fig. [6] but now we focus on the compact five planet case with the initial configuration 
3:2, 3:2, 4:3, 4:3 and a planetesimal disc of 35M^. The outcome is similar to, but slightly better 
than, the resonant Nice case with four planets: we have 6 cases where the giant planets are in the 
right positions and have their current obliquities, but in only 2 of these Arojs circulates. There are 
fewer cases where Jupiter is tilted more than Saturn: ej > es 37% of the time for all simulations 
that produce four planets, and ej > es 33% of the time for the cases when the planets are close to 
their current orbits. But (ej/es) = 1.44 ± 2.12 for all final cases with four planets, and 1.45 ± 2.02 
for when the final planets are close to their cnrrent orbits, again indicating that there is a longer 
tail at large values of ej/es. 

Finally, we turn to the loose five planet system, whose initial confignration is 3:2 3:2 2:1 3:2 
and the ontermost ice giant resides at 22.2 AU. The mass of the planetesimal disc is 20M^. We 
plot the final obliquities of Jupiter and Saturn in Fig. [8l The main difference with the same figures 
for the resonant Nice case (Fig. [6]) and the compact five planet confignration (Fig. [7j) is an absence 
of high es a-t low ej. Only one simnlation with fonr planets too far from their correct positions at 
the end has final es > 15° and ej < 5°. The fraction of simulations with ej > es is 45% for all 
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Fig. 8.— Same as Fig. [H but for the A^-body simulations of the loose (3:2 3:2 2:1 3:2) five planet 
model. 

runs that have four planets at the end and 61% for runs where the planets are near their current 
orbits. We also find that (ej/es) = 3.98 ± 3.45 and 4.39 ± 7.24 for both samples, respectively, so 
that this wider initial configuration of the planets leaves Jupiter’s obliquity higher compared to 
Saturn’s than the compact one. 


5.2.2. Evolution examples 

We now illustrate the mechanisms behind the tilting of Jupiter and Saturn through some 
figures. As was discussed earlier, tilting Saturn requires that as resonates with ss, or passes 
through this resonance. For Jupiter, on the other hand, there are two ways to tilt it: passage 
through resonance with sg or sy. Figure [9] shows the evolution of the system that tilts both Jupiter 
and Saturn to their current obliquities in a simulation of the compact five planet model. The 
top-left panel shows the perihelion distance, semi-major axis and aphelion distance of all planets 
as a function of time. To guide the eye each planet is identified by a different colour. The top-right 
panel depicts the period ratio between Jupiter and Saturn. The bottom-right panel depicts ej 
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Fig. 9.— Evolution of a compact five planet case, where Saturn gets tilted and Jupiter does not. 
The top-left panel depicts the semi-major axis, perihelion and aphelion of the planets versus time. 
The letters and colour coding guide the eye. The top-right panel depicts the period ratio between 
Jupiter and Saturn versus time. The grey band is the region between 2.1 and 2.3 that needs to 
be crossed quickly. The dashed line shows the current value. The lower-left panel shows the spin 
precession constants of the gas giants and the vertical eigenfrequencies sj and sg as a function of 
time. Last, the bottom-right panel shows the obliquities of the gas giants with time. 
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Fig. 10.— Same as Fig. [U but for a case where Jupiter gets tilted but Saturn does not. 
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(black) and es (blue) and the bottom-left panel shows the evolution of aj (magenta), as (black), 
|s 7 | (red) and jssj (blue) as a function of time. When starting with three ice giants there is an 
extra vertical eigenfrequency that disappears when the third ice giant is ejected. The frequencies 
decrease with increasing semi-major axis and this property allowed us to keep track of which 
frequency corresponded to which planet and ultimately identify sy and sg belonging to the two 
surviving ice giants. 

After approximately 3 Myr of evolution, the system becomes unstable. The innermost ice giant 
becomes Neptune, and the outermost ice giant evolves little and ends up being Uranus. The period 
ratio between the gas giants increases rapidly when Saturn scatters the lost ice giant and then 
evolves more slowly. The two remaining ice giants eventually settle into stable orbits at f ~ 7 Myr 
and migrate towards their current orbits without undergoing encounters. The frequency Isgl crosses 
aj shortly after the system goes unstable but the crossing occurs during the jump and is too quick 
to allow capture in resonance. It then steadily approaches as and after about 100 Myr it is close 
enough for Saturn to be caught in resonance. Saturn’s obliquity increases and it stays in resonance 
until the end of the simulation. The resonant state is marked by the long-period oscillations in 
the obliquity of Saturn. At the same time Isyj approaches aj but does not come close enough for 
capture in resonance because Uranus stays fairly close to the Sun. 

From Section [2] we know that a slow migration of Neptune is needed to tilt Saturn. From 
Fourier analysis we obtain for the above simulation Igg = 0.065° and as = 0.693" yr“^, where Igs 
is the forced inclination in Saturn associated with sg- Thus Ts^min ~ 10 ® yr, according to Figure [T ]6 
and equation ([2|). In the simulation Neptune migrates outwards by 0.4 AU between t = 80 Myr and 
t = 350 Myr, and then essentially stops. This suggests that a Aa/At ~ 2 X 10 ® AU yr The 
time scale for migration is then tmig = a/a ^ 10® yr and thus Saturn gets captured in resonance. 

Unfortunately we also have many cases where Jupiter gets tilted significantly more than Saturn. 
In Fig. [10] we plot the evolution of a simulation of the compact five planet model where Jupiter 
gets tilted but Saturn does not. The panels are the same as for Fig. [H The outermost ice giant is 
ejected when the system becomes unstable at f ~ 9 Myr. From the bottom-left panel we see that 
aj gets caught in resonance with sj. Thus it appears that cases of high obliquity of Jupiter are 
often caused by capture or passage through resonance with sy rather than sg because the latter 
occurs early in the violent phase when the planets migrate fast. Applying Figure [T]& to Jupiter, we 
have I 57 = 0.072° and aj = 3.24" yr“^ and thus Ts^min ~ 200 Myr. Uranus migrates outward by 
0.3 AU from f = 50 Myr until t = 200 Myr so that tmig 3> 200 Myr. Thus Jupiter gets caught in 
the resonance. 


5.2.3. Spin-orbit resonance erossing 

Figure [11] shows the cumulative distribution of Jupiter’s final obliquity, compared with the 
cumulative distribution of Jupiter’s obliquity before any |s 7 |/aj crossing (if one occurs), for the 
simulations with four planets in the right positions at the end (blue dots in Figs. [BHS]). Except for 
one simulation in the compact five-planet model, all simulations have ej < 10 ° before any |s 7 |/aj 
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Fig. 11.— Cumulative distribution for Jupiter’s obliquity to be less than ej for the A^-body simu¬ 
lations with four planets in the right positions at the end. The resonant Nice, compact five-planet, 
and loose five-planet models are shown in red, black, and blue, respectively. The filled symbols 
show the cumulative distribution of the final obliquity, while the open symbols show the cumulative 
distribution of the obliquity before any |s 7 |/aj crossing (if one occurs). 

crossing occurs, but at least 18% of the simulations have final ej > 10°. This confirms that the 
cases of high Jupiter obliquity are often caused by capture or passage through resonance with sj 
instead of sg- 

To understand the differences in the distributions of the final obliquities of Jupiter and Saturn 
between the three models, we compare in Figures [T2Hin the final values of |s 7 |, jssl, aj, and as for 
the A^-body simulations with four planets in the right positions at the end. In each figure, the filled 
circles and error bars show the mean values and three standard deviations (3a) of the final values 
of |s 7 | versus au and [sgl versus on, with 57 and sg from the Laplace-Lagrange secular theory. The 
red line shows the final aj, which is fixed after rescaling the final aj to 5.20 AU, and the cyan band 
shows the ±3 (t range of the final values of as. The vertical magenta dashed lines indicate the ±10% 
range of the current au and on- 

From Figure [m for the resonant Nice model, we see that |s 7 | is about 0.7a above aj, which 
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Fig. 12.— Final values of |s 7 | versus au and jssl versus on for the Al-body simulations with 
four planets in the right positions at the end of the resonant Nice model (blue dots in Fig. [6]), 
with the filled circles and error bars showing the mean values and three standard deviations (3 cj), 
respectively. The red line shows aj, and the cyan band shows the ±3cr range of the final values of 
as- The vertical magenta dashed lines indicate the ±10% range of the current au and on- 


means that jsyl crosses aj in about 24% of the simulations with four planets in the right positions, 
if we assume a normal distribution. This agrees with the ~ 24% of simulations with final ej > 10° 
(see Fig. [m, showing again that the cases where Jupiter has a high obliquity are most likely 
caused by a resonance crossing with sy instead of sg. For the compact five planet model (Fig. [T3]) . 
|s 7 | is further (~ la) above ctj, and there are fewer (~ 20%) simulations with ej > 10°. On the 
other hand, although more than half of the loose five planet simulations have final |s 7 | < aj (see 
Fig. [IT)) , only ~ 18% of the simulations have ej > 10°. This suggests that, in this model, the 
|s 7 |/q!j resonance crossing often occurs when Uranus is migrating too quickly, which is consistent 
with the fact that Uranus is on average farther from the Sun than in the other two models. 

Unfortunately, this problem of resonance crossing occurring too quickly also prevents the tilting 
of Saturn in many cases. For the resonant Nice model (Fig. fT2]) . more than half of the simulations 
have IssI < <aS) but only ~ 8% of the simulations have ss > 10°) indicating that the |s8|/<as 
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Fig. 13.— Same as Fig. [T^ but for the A^-body simulations with four planets in the right positions 
at the end of the compact five planet model (blue dots in Fig. [7]). 


resonance crossing often occurs when Neptune is migrating too quickly. Although the compact 
five planet model has a similar mean value of on, it has a narrower range in un and more overlap 
between jssl and as (see Fig. n, which results in a higher fraction (~ 19%) of simulations with 
es > 10°. Finally, in the loose five planet model (Fig. fT^ . Neptune is on average too far from 
the Sun, and almost all of the simulations have jssj < as, but only ~ 3% of the simulations have 
es > 10°. 


5.2.4- Summary 

We summarize the outcome of our N-body simulations in Table [TJ The column is the 
probability of the simulation having four giant planets after 500 Myr. The column fa is the 
probability that the system has four planets with each planet within 10% of its current orbit. The 
next column, f^, is the probability of having four planets close to their current orbits and Arojs 
circulating. The last column, f^, is the probability that the system also has ej < 5° and es > 15°. 
For the resonant Nice and compact five planet models, we list two values for /g, with the lower 
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Fig. 14.— Same as Fig. [T^ but for the A^-body simulations with four planets in the right positions 
at the end of the loose five planet model (blue dots in Fig. [8]). 


value requiring the circulation of Arojg and the higher value removing this restriction on the state 
of Atujg. For the loose five planet model, none of the simulations with four planets close to their 
curr ent orbits has e.^ < 5° a n d eg > 15°, and there is only an upper limit on f^. 


Nesvornv fc Morbidellil ([2013) have examined the statistics of orbital properties for a large 


number of models, each with only 30-100 A-body simulations. In this paper, we have examined 
three representative models with 1280 simulations each. Our resu lts on the orbital statis t ics ar e 
much more accurate, and they are consistent with those found bv iNesvornv &: Morbidellil (|2012l b 
In particular, the loose five planet model is the most successful in reproducing the positions of 
the giant planets and the circulation of Atujs. However, it is clear from Table [T] that the loose 
five planet model fails to reproduce the obliquities of Jupiter and Saturn (with the probability 
fe < 0.08%) with the setup we have chosen. Thus the obliquities of Jupiter and Saturn provide a 
strong constraint that can alter the relative merits of different models and initial conditions. 

In summary, it appears that both the resonant Nice and compact five planet models are able to 
reproduce simultaneously the current orbital architecture of the giant planets and the obliquities of 
the gas giants, with the compact five planet model faring slightly better. However, the probability 
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Table 1: Statistics of A^-body Simulation Results 


Model 

-^disc 

Initial on 

k 

fa 

fuj 

fe 

3:2 3:2 4:3 

50 M® 

11.6 AU 

35% 

6.6% 

1.0% 

0.08-0.23% 

3:2 3:2 4:3 4:3 

35 M© 

14.2 AU 

23% 

6.7% 

1.9% 

0.16-0.47% 

3:2 3:2 2:1 3:2 

20 

22.2 AU 

27% 

10% 

3.3% 

<0.08% 


is less than 0.5%, even if we lift the restriction on the state of Arojs- The loose five planet model 
fails with a probability of less than 0.08%. 


6. Discussion 


Even though the results presented in the previous section are extensive, they are by no means 
complete. In this section we discuss a few outstanding issues. 


First, we did not impose any condition on the migration speed of the giant planets. iNesvornv fc Morbidelli 


(j2012l ) imposed having the Saturn-Jupiter period ratio evolve from below 2.1 to beyond 2.3 in 
shorter than 1 Myr (their criterion D). We verified that all nine cases in our A-body simulations 
that reproduce the current obliquities of the gas giant have some sort of jump in the period ratio, 
although only one case from the compact five planet model had a jump in the right range, and 
is depicted in Fig. [9l This would suggest that the compact five planet case is better overall at 
matching all constraints than the resonant Nice model. 

Second, we only tested a single four planet and two five planet configurations because we opted 
to bracket the two end points of the initial spacing of the giant planets to determine whether the 
final outcomes favour one case over the other. The ai iswer is clearly yes. W e have had to run 
many more simulations per set of initial conditions than INesvornv fc Morbidellil ([2013) because we 
are trying to satisfy an additional constraint, namely reproducing the obliquities of the gas giants. 
Given that we found at most a few cases out of 1280 that match all constraints from one set of 
initial conditions, it becomes evident why we ran as many simulations as we did for each set of 
initial conditions, and why it is infeasible to test a larger sample of initial conditions. 

Third, we generally did not change the mass of the planetesimal disc when running the simu¬ 
lations for one configuration. One could argue that if we increased the disc mass from 35 to 
50 in the compact five planet case we may get a higher probability matching all constraints. 
However, increasing the disc mass from 35 to 50 would most likely damp the eccentric- 
ities of Jupiter and Saturn to o much and make them inconsistent with their current values - see 
Nesvornv fc Morbidellil ([20121 ) for a more detailed discussion. We ran the same compact five planet 
model with a 50 Mq disc. Even though the probabilities /4 to show no improvement and are 
nearly identical to those of the compact five planet model with the 35 disc, we found that the 
eccentricities of Jupiter and Saturn were inconsistent with their current values (> 2cr result). Thus, 
increasing the disc mass does not help our case and we do not report the results here. In addition, 
our current disc mass of 35 typically increases the period ratio of Jupiter and Saturn by 0.2 
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after the jump, which iNesvornv &: Morbidellil (|2012l j state may indicate the disc is too massive. 

Anoth er issue that we have not consider ed is the uncertainties in the spin precession constant 
aj and ag. IVokrouhlickv &: Nesvornvl (j2015l i have published an analysis based on one-planet sim¬ 
ulations similar to those in Section [2] (i.e., integration of spin evolution for a planet whose orbital 
nodal precession rate |s| is decreased on timescale r^). Their analysis differs from ours in assuming 
that the hnal values of sy and sg are hxed at the observed Solar System ones, but they considered a 
range of values for both aj and ag. For the Isgl/aj resonance crossing, they hnd constraints on hs 
and Ts for the resonance crossing to be too fast for capture and Jupiter’s obliquity to be < 3° after 
this crossing (which agrees with the open symbols in Fig. ITT]) . Similarly, for the final approach of 
sy to its current value, they find constraints on og and that yield Jupiter’s final obliquity to be 
< 3°. For the final approach of sg to its current value, they considered different values of ag and Tg 
and find that capture into the |sg|/Q;g resonance requires large and low gg (inc luding the value 
of ag adopted in the present paper). However, IVokrouhlickv fc Nesvornvl (j2015l l did not assess 
the probability that su c h con ditions on the migration rates are satisfied in the various models of 


Nesvornv &: Morbidellil (j2012l l. nor did they assess the probability of these models giving the right 


values of S 7 and sg, which we have done here for three specihc conhgurations. Thus, while they 
suggest that the mechanism behind tilting Saturn is the resonance crossing between sg and as, 
we show here how well this mechanism performs with full V-body simulations in several specihc 
settings. 

A last point that requires elaboration is whether the obliquities of the gas giants can be used 
to constrain the migration of the giant planets at all. At issue is the following. The current solar 
system has £j ~ 3° and eg ~ 27° and for both planets aj cosej ~ Isyj and ag cos eg ~ |sg|. With a 
10 % margin in the hnal semi-major axes of the planet one can compute the range of expected hnal 
obliquities of the gas giants assuming that during a resonance crossing there is perfect capture i.e. 
cosej = Is^/ojI if Is^/ckjI < 1, and cosej = 1 otherwise. Here primed quantities are the values 
at the end of the simulation. Doing a series expansion of the Laplace-Lagrange theory around the 
current semi-major axes of the ice giants we hnd that e'j G (0,44°) and eg € (0,53°). With such a 
large range of hnal obliquities over a small range of hnal semi-major axes of the giant planets, one 
might argue that the hnal obliquities of the gas giants cannot be used to constrain the migration. 

However, this argument requires one assumption that does not stand up to scrutiny. Indeed, 
if there was perfect capture then in a plot of \s/a\ vs e the data points should be clustered near 
e = arccos(|s/a|) when \s/a\ < 1 and near e = 0 otherwise. Figures [T5][T7] depict the hnal obliquities 
of the gas giants as a function of | s/a|. It is clear that this trend is only partly visible. We see 
some clustering of high sj near |s 7 /aj| ~ 1 - and in some cases near |s 7 /aj| > 1 because of the 
inherent uncertainty in the Laplace-Lagrange model. For Jupiter the trend is more obvious than 
it is for Saturn, and the trend for Saturn is almost absent in the loose hve planet model. This 
indicates that when a crossing does occur, capture often does not. In other words, the assumption 
that during a crossing there is perfect capture is clearly untrue. 

So what are we to make of this? The answer is that for Saturn’s obliquity to be where it is, two 
things need to be satished simultaneously: a) Neptune’s migration needs to be slow enough so that 
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Fig. 15.— The final obliquity of Jupiter (left) and Saturn (right) vs Isy/ajj (left) and Iss/asI 
(right) for the resonant Nice model. The dashed line shows e = arccos(|s/Q:|). 

r X i > 30 Myr deg, and b) the spacing of the planets needs to be such that Isg/asI ~ 1. We showed 
in Figs. [T^14l that the final semi-major axes of Uranus and Neptune (and indirectly, also the one of 
Saturn) are not the same for each set of simulations. Consequently, the final values of sj and sg are 
not the same either, despite the requirement of the final orbits being within 10% of the current ones. 
As an example, for the resonant Nice model, the average and standard deviation of the final semi¬ 
major axes of Saturn, Uranus and Neptune in the cases where the planets are near their current 
orbits, are (as) = 8.98T0.352 AU, (ou) = 18.69±0.99 AU and (on) = 30.14±1.70 AU. In contrast, 
for the compact five-planet model the values are (as) = 9.33 ± 0.451 AU, (au) = 18.69 ± 0.956 AU 
and (un) = 29.74 ± 1.08 AU. Even though these are all within 10% of the current values, the 
current spacing of the giant planets is better matched for the compact five planet case than for 
the resonant Nice case. It is the difference in the final spacing, coupled with the tail of Neptune’s 
migration, that determines the obliquities of the gas giants and whether one model matches the 
observed constraints better than another. In light of this, we conclude that the obliquities of the 
gas giants can indeed be used to constrain the evolution and initial conditions of the giant planets. 
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Fig. 16.— The final obliquity of Jupiter (left) and Saturn (right) vs Isy/ajj (left) and Iss/asI 
(right) for the compact 5 planet model. The dashed line shows e = arccos(|s/a|). 


7. Summary and conclusions 

We have performed secular and Wbody simulations tracking planetary obliquity to study 
constraints on migration histories of the outer Solar System from both spin and orbit, covering the 
smooth migration scenario, the resonant Nice model and two five planet scenarios. The obliquities 
of Jupiter and Saturn provide a strong constraint on the models. The secular results leave us in 
the fo llowing situation: in the very simplest case, that of smooth migration, the lWard fc Hamilton 
(j2004l ) secular spin-orbit resonance mechanism for tilting Saturn’s spin axis appears to work nicely 
if the product of the migration time scale and the orbital inclinations is sufficiently large {t x i > 
30 Myrdeg). On the other hand, the resonant Nice model, which is preferable on many grounds 
(such as explaining the orbital properties of the giant planets, terrestrial planets, and main belt 
asteroids), is more problematic, because the secular eigenfrequency sg is initially above the spin 
precession frequency of Jupiter, uj. The fundamental problem is that there are incompatible 
constraints. We need to migrate slowly (r > 8Myr) at typical inclinations to consistently tilt 
Saturn’s spin axis by > 15°, but for those same inclinations we need to migrate quickly (say, 
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Fig. 17.— The final obliquity of Jupiter (left) and Saturn (right) vs |s 7 /aj| (left) and |s8/«s| 
(right) for the loose 5 planet model. The dashed line shows e = arccos(|s/a|). 

r < 2Myr) to avoid tilting Jupiter by more than the observed 3°. We find that on average ej ~ es 
in our secular simulations of the resonant Nice model. 

At the same time, the Wbody simulations appear to tell a different story. The resonant Nice 
model is able to reproduce the current orbital architecture of the giant planets and the obliquities 
of Jupiter and Saturn, but only with a small probability. The compact five planet model fares 
slightly better, but the loose five planet model has great difficulty reproducing the obliquities, even 
though it is the most successful in reproducing the orbital properties. However, it is possible that 
a slight change in the initial conditions or the outer edge of the planetesimal disc could improve 
the outcome. 

Ultimately the following needs to happen; (1) There needs to be fast migration during the 
encounter phase to avoid tilting Jupiter through resonance passage with sg- (2) Then a late, slow 
migration of Neptune to its current location could complete the task by tilting Saturn through the 
resonance with sg- (3) At the same time Uranus must stay close enough to the Sun to avoid tilting 
Jupiter through the resonance with sj. Condition (1) is almost always satisfied in the Wbody 
simulations. But the chances of satisfying conditions (2) and (3) are limited, as the jssl/as reso- 
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nance crossing often occurs when Neptune is migrating too fast, especially in the loose five-planet 
model, while Jupiter sometimes crosses the |s 7 |/q;j resonance. Thus the main obstacle encountered 
to reproduce the obliquities of both Jupiter and Saturn is the final orbital spacing of the giant 
planets, coupled with the tail of Neptune’s migration. Both the resonant Nice and compact hve 
planet models appear to be able to overcome the obstacle, but with low probability. 
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